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— ^ Abstract 

Q We analyze in detail a one- dimensional stochastically driven running sandpile. 

O, The dynamics shows three different phases, depending on the on-site relaxation 



rate and stochastic driving rate. Two phases are characterized by the presence of 
Cd travelling waves. The third shows algebraic relaxation. 
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1 Introduction 

The investigation of a simple automaton with nontrivial time correlations by Bak, 
Tang and Wiesenfeld [1] has triggered a number of studies on a variety of similar 
models. Whereas the initial model was deterministically driven at a vanishing 
rate (every perturbation was allowed to relax to a stationary state before a new 
perturbation was added), later models (for instance [3, 4]) allowed for stochastic 
driving at a finite rate. The case of infinitely extended lattices has been analyzed 
using methods from renormalization and mode coupling theory ([4]). For finite 
lattices, multifractal scaling has been used to analyze finite size effects ([2]). 

We here analyze for a specific automaton both the deterministic dynamics and 
the stochastic relaxation in some detail. Although there is no system which is ex- 
actly described by this automaton, it helps to have in mind some physical processes 
with qualitatively similar behaviour. As in the Bak, Tang, Wiesenfeld [1] model, 
one may interpret our automaton as an idealized one- dimensional sand pile. Be- 
cause of the boundary conditions, it is more like a sand pile near a wall. Sand is 
added randomly at a site, which topples if the slope exceeds a certain threshold. 
Alternatively, one may think of a chain of water reservoirs, where water is added 
randomly and where the reservoirs spill over, if a certain water level is exceeded. 
We assume that the amount of mass transported per time step in the case of a 
supercritical slope consists of a constant part plus an additional slope dependant 
part. One expects and finds in numerical simulations that, after some transients, 
the automaton relaxes to a statistically stationary state with nontrivial dynamics. 

Analysis of this stationary state reveals three different kinds of phases, depend- 
ing on the driving rate and the relaxation rate. In the Bak, Tang, Wiesenfeld [1] 
case of infinitesimal driving rate the dynamics is rather trivial and consists of waves 
which pass through the system with rather little interaction. As the driving rate is 
increased, there is first a regime dominated by a different kind of waves, and then 
a regime with algebraic relaxation because of interaction between relaxation and 
driving. The boundaries between the different regimes can be determined and are 
in agreement with numerical simulations. 

In section 2 we define the model and its parameters and present simulation 
results for the slope profile of the statistically stationary state, for the dynamics of 
perturbations and for the temporal correlations at a fixed site. 

In section 3 we discuss some idealized solutions to the dynamical rules which 
may be thought of as waves and wave packets. In some regions of the lattice they 
exist as undamped running waves whereas in others they are exponentially damped. 
A transport argument allows to explain quantitatively the structure of the slope 
profiles. 

In section 4 we use methods from stochastic calculus to obtain a formally exact 



solution for the relaxation dynamics. This allows us to discuss not only the low 
frequency behaviour which is not accessible in direct numerical calculations, but 
also to estimate the transition between the various scaling regimes in the power 
spectra. 

In section 5 we compare our results with previous ones and comment on their 
relevance to other, in particular higher dimensional automata. 

2 Definitions and simulation results 
2.1 The model 

As described above, one may think of our model as a chain of water reservoirs 
or as a one- dimensional sand pile. This picture suggests to use the water level at 
each reservoir or the height of the sandpile at each site as variables. However, the 
dynamics is controlled not so much by the absolute level but rather by the amount 
of mass transported, so that the slopes, i.e. differences between neighboring sites, 
are better suited variables for these automata. 

The state of the system may thus be described by a vector h = (ho, ..., hjM_i) G 
i?^ for the heights or by a vector a = ((Tq, ..., (Tat-i) G i?^ for the slopes. The 
interpretation of these vectors can be seen in Fig. 1. The transform 

h = (ho, ...,hN_i) I — > a(h) = (do, ...,aN_i) = (/iq - /^i, ..., /iAr-2- /^w-i, /iw-i) (1) 

maps between the two description pictures. Henceforth we will use the slope picture. 
The dynamics consists of two steps: A relaxation rule and a perturbation rule (the 
addition of mass). The state a at time t evolves into the state a' at time t + 1 
according to these steps: 

(i) Relaxation rule 

H a slope (T,- is too steep it will relax by a certain amount r,- of mass being 
transferred to the right neighbouring site, r,- is given by 



r,: := < 



1 + (ai - (Tc)/P if (Ti > (Tc 

else 



, i = 0,...,A- 1. 



One arrives at an intermediate state ct P: 
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The relaxation rule contains two parameters, a critical slope a^ and a relax- 
ation rate P. 

(ii) Perturbation rule 

A site ig of the lattice is chosen randomly and a mass s G -R"*" is added to that 
site. For the first site, ig = 0, we have 



0,...,iV-l (4a) 



and for all other sites, ig 7^ 0, 

tmp . , . . . 

o[:=\ a^tmp^^ ^ ^^^^ , ^ = 0,...,iV-l . (4b) 

tmp . . -, 

o^ ^ - s , i = io - 1 

Here another parameter s > characterizing the driving rate enters. The 
perturbation rule is illustrated in Fig. la and lb. 

We expect the model to evolve from an arbitrary initial state into a statistically 
stationary state. In order to recognize the time when the model reaches this state 
we averaged the slope profiles in blocks of lOOiV time steps and calculated the 
relative differences in the profiles between consecutive blocks. If this relative error 
became less then a prescribed value we assumed the model to have reached the 
stationary state. To be on the safe side averages were computed after an additional 
iteration of 10 times the number of steps needed so far. The value of the critical 
slope (Tc only affects the average slope, i.e. the height of the stationary state in 
the slope picture, but it has no effect on the structure of the stationary state or 
dynamical properties. The essential parameters are therefore P and s; for the trivial 
parameter a^ we chose CTc = 4 throughout the paper. In the next two subsections 
we illustrate the properties of the stationary state. 

2.2 Slope profiles in the stationary state 

The model relaxes to different stationary states depending on the two parameters 
P and s. Fig. 2, a division of the P-s-plane which separates automatons with 
structurally different stationary states, directs attention to Fig. 3-7 in which time 
averaged slope profiles along with some snapshots out of certain lattice parts are 
shown. Common to these figures is the choice for the lattice size N = 200 and the 
critical slope a^ (the critical slope is represented by dashed lines in the figures). We 
have omitted plotting the average slope at the last lattice site in each profile since 
it is affected by a trivial boundary effect. 



One can infer from these figures the apparent importance of wave packets which 
are visible in all of the snapshots but frame 3 in Fig. 7. We can distinguish packets 
with averages below and above the critical slope a^. Following the time evolution 
of these profiles one notes that the wave packets are stable and propagate one site 
per time step. If a wave packet is perturbed by the addition of mass the resulting 
propagating profile is still a wave packet but it can contain a different number of 
periods than before. Wave packets arise from perturbations in parts of the lattice 
with constant slope. The exceptional case mentioned above (part 3 of the profile for 
P = 6 , s = 2) exhibits another kind of behaviour: The addition of mass results in 
a perturbation of the straight line which will be damped during the time evolution. 

We refer to different parts of the stationary state as phases characterized by the 
excitations they carry: Subcritical wave packets in subcritical regions, supercritical 
wave packets in supercritical regions and damped perturbations in linearly increas- 
ing lattice parts. The division of the P-s-plane in Fig. 2 separates automatons 
with different numbers of phases: The ones in region I contain only the phase with 
subcritical waves, the ones in region II contain the first two phases and the ones in 
region III all three phases. 

2.3 Temporal correlations 

Temporal correlations in regions of linearly increasing average slope (phase 3) show 
algebraic decay as demonstrated in Fig. 8 for the parameters P = 10 , s = 2. 
Fig. 8a shows the autocovariance function of the slope at a fixed position in this 
phase. The spectral density function plotted in Fig. 8b consists of two clearly 
distinguishable regions of algebraic behaviour with different slopes. Increasing the 
value of P leads to higher scaling exponents. 

In the next two sections we attempt to explain the simulation results. 

3 Wave packets and transport 

Wave packet like structures are characteristic excitations in the first two phases. 
As will be shown in subsection 3.1 the dynamical rule of the model supports certain 
exact running waves and wave packet like solutions which are undamped. In con- 
trast, excitations in the third phase are damped as will be proved in subsection 3.2 
thus explaining the absence of wave packets in linearly increasing regions of the 
lattice. A transport argument given in subsection 3.3 based on these excitations 
finally explains the partition of the parameter plane. 



3.1 Idealized waves and wave packets 



We consider in this subsection a simplified system with an infinite number of sites, 
i.e. the state space of the system is the set of series ((Jn)nez- The dynamic reduces 
to the relaxation rule 



Ti := < 



1 







, else 



a'- = Oi + ri_i - 2ri + r^+i , for alH G ^ . 
An idealized infinitely extended wave is a state 



(5) 



(6) 



V 



c + (i , n even 
c — d , n odd 



\ 



(7) 



/ nGZ 

characterized by two parameters c £ R and rf > 0, whose evolution in time is simply 
a shift of one place per time step. A subcritical wave packet is one with an average 
below criticality, i.e. c < a^ and c + d > a^; a supercritical wave packet is one with 
an average above criticality, i.e. c > a^ and c — d < a^. The pairs (c, d) G i?x]0, oo[ 
for which an idealized wave obeys the dynamical rule in eq. (6) obviously depend 
on the parameter P through eq. (5). One finds for 

• subcritical waves 



P > 1 : T^gub = {(c, d) \ a, - 1< c < a, , 

^ = 1 : l^sub = {(c, d) \c = a,-l, d>l} 

0<P< 1 : W^^^ = {{c,d) |c<c7, -1, 

• and for supercritical waves 



P > 2 

1< P < 2 
< P < 1 



VFsnner = {(c, d) \ a,. < C < a,. + -^ 



'super 



rf=l + ^±iff-} 



l^super = {(c, d) \a,<c,d=l+ ^±^} 



W, 



super 



(9) 



In the subcritical case it is easy to see that besides these infinite structures wave 
packets with finite extent are also solutions of the dynamical rules. A wave packet 
of finite extent has a certain number of periods (c — d,c + d) or (c + d,c — d), 
bracketed by infinite sequences of c, for instance 



a = {. . .,c,{c- d,c + d)„,c,. . .) 



:io) 



This structure (10) propagates one site to the right per time step. 

Finite wave packets are results of perturbations in constant parts of the lattice: 
It is for instance possible to observe in the time development of the state 



o- = (..., c, c - s, c + s, c, ...) , c< a,. , c + s > a,. , 



111 



obtained by adding mass to a certain site, the formation of a finite wave packet 
where the number of periods depends on the choice of the parameters P, s and the 
mean value c. 

In the supercritical case things are a little more complicated: There exists no 
state with only a finite number of components deviating from a constant c which 
propagates without changing its form. But it is possible to construct propagating 
structures which are similar to ideal wave packets of finite extent. Consider for 
instance a state with the properties 

Vi < : (Tj- = c > (Tc 

0-1 < (Tc 

0-2 > (Jc 

0-3 < (Jc 

Vi > 4 : (Tj- > (Tc 

A necessary and sufficient condition for this state to propagate backwards one site 
per time step without changing its form requires special values for the components 
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13) 



The possibility of wave packets with an arbitrary number of periods in the sub- 
and the supercritical case is of great importance for the transport argument below. 



3.2 Damping of perturbations 

We consider a part of the lattice with always supercritical slope. Then the relax- 
ation rule simplifies to 



Ti 



1 



(Jr 



P 



a„ a,- 

1 £ + ^ 

P P 



, 1 



P 



P 



I 1 



+ 1 



(14) 



(15) 



We are particularly interested in states which depend linearly on the index i because 
the time averaged slope profile in part 3 of the lattice shown in Fig. 5 is of that form. 
One property of the relaxation rule (15) is that it leaves straight lines ai = c + ib 
invariant. Localized perturbations of these straight lines expand only one site per 
time step and thus remain localized. For the evolution of such perturbations it 
then suffices to study the dynamics on a finite sublattice of size N . Because of 
linearity of the dynamics described by (15) we may separate the evolution of the 
perturbation e(n) G R^ (n is the time step) from an underlying time independent 
straight line and write down a vector equation 



s{n + 1) = Ms{n) = M"+^e*(0) 



(16) 



where 



M 



P-2 
P 


1 
p 





1 
P 














V 
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1 P-'2 



i, , 



(N X N - matrix) 



(17) 



/ 



is the relaxation matrix. Since M is symmetric it may be diagonalized. Its eigen- 
values are 

P -9 9. h.TT 

18 



l^k 



P-2 2 kir 

^r cos ^^ r , k = l,...,N 



P P N + 1 

with the corresponding eigenvectors 



iW 



-1 



\n + l 



nkiT 



■ sm ■ 



, n=l,...,iV . 



(19) 



\l N + 1 N + 1 
After that diagonalization the derivation of a formal expression for e(n) leads to 



N 






f20) 
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Since for P > 2 all eigenvalues have an absolute value less than one, any perturba- 
tion will be damped. For the model this means that in parts of the lattice which 
are always supercritical all deviations from a straight line tend to vanish during 
time evolution. 

3.3 Transport argument 

The mean transport rate at site i of the lattice must be (i + 1)-^ since the average 
mass added per time step at a fixed site is s/N . In order to handle this fiux time 
evolutions with the right amount of excess slope are necessary. 

Consider first the case P > 2. Subcritical wave packets separated by regions of 
slope less than a^ cover a transport rate range of 0,|-^3y since the transport in 
regions of undercritical slope vanishes, while in a region containing a wave packet 
with average c = a^ it is maximal: Even sites contribute r,- = -^^j ^^^ ^^^ sites 
nothing, resulting in an average transport rate |-^rY' From these considerations we 
conclude that for s < ^jrri the transport can be achieved through subcritical wave 
packets. An important point is the variability in the length of the wave packets and 
in the subcritical regions between them which leads to a variety of time evolutions 
covering the whole transport rate range 0, |"^rY ■ ^^^ ^^^ °^ these evolutions have 
the same time averaged slope profile: A horizontal subcritical line. This explains 
the plateaus seen in Fig. 3. 

For s > ^jTTi time evolutions with transport rates greater than the rates of 
subcritical waves must enter the dynamics. Larger rates can be achieved by su- 
percritical wave packets separated by parts of constant supercritical slope. The 
minimum rate is attained by a packet that extends over the entire lattice region 
considered since now only every second site is supercritical. Its transport rate is 
exactly 2'p^- "^^^ state with the maximum rate has a region with constant su- 
percritical slope. The maximum slope in such a region is a^ + j^ because we still 
want to retain the possibility of supercritical wave packets. Thus transport rates in 
IpTT' p^ ^^^ ^^ achieved by states consisting of supercritical waves separated 
)y regions of constant supercritical slope. 

Therefore we expect in the case P > 2 and ^jtZi < •§ < ^^ a stationary 
state consisting of two regions: A constant subcritical region proceeding from the 
left boundary and a constant supercritical region extending to the right boundary. 
The transition occur approximately at the site at which the needed transport rate 
exceeds the maximum rate that can be achieved by subcritical waves. From this 
one estimates that the transition between region 1 and 2 should occur near 

1 P A 

rai^2 = 1 • (21) 

2P- 1 s ^ ^ 

These considerations are also in agreement with the simulation results as can be 
seen in Fig. 4. 

9 



At parameter values s > -j^ supercritical waves do not suffice to provide the 

necessary transport. This implies the existence of a third region in the stationary 

state which is different from the first two regions. The excess slope is so large that 

no wave packets are possible due to the damping of perturbations in this part. 

Therefore, there can be no variability in transport rates with time evolutions with 

the same time averaged horizontal profile. We expect in this case an increasing 

slope in order to adjust the transport. The transition point between second and 

third region can be estimated by an argument analogous to the one utilized above 

to be 

P -IN 

"-» = p-2T - ' ■ ('"^' 

Furthermore, one can show that the slope of the profile in the third region must 
be Pj^ and the averaged value at site N — 1 must be < (Jn-i >= Pis — 1) + (Jc 
(because of fiux continuity). All these predictions are consistent with simulation 
results (cf. Fig. 5). 

For the case P g]1, 2[ analogous considerations apply. Subcritical waves cover a 
transport rate range 0,|-^3y while every rate larger than ^jrri can be achieved 
by supercritical wave packets. For s < ^jrri we expect the stationary state to be a 
subcritical line which is confirmed by the results depicted in Fig. 4. For s > ^jrri 
a region with time averaged slope greater than a^ is necessary. In this region 
supercritical packets are characteristic of the dynamics. An estimation of the site 
at which the profile crosses the critical value is again 

1 P N 

ni^2 = 1 • (23) 

2P- 1 s ^ ^ 

Since for these values of the transport parameter P there is no limitation for the 
mean value c all transport with high rate can be done by supercritical wave packets. 
But in order to adjust the transport rate to increasing (with site index i) demands 
we expect the mean value c to increase also with the site index i. These estimates 
are in good agreement with simulations (cf. Fig. 7). 

The arguments given above finally lead to a partition of the P-s-plane in three 
distinct parts separated by the curves s = ^jtti ^^^ s = ^^ as shown in Fig. 2. 



4 Calculation of temporal correlations 

In this section we concentrate on the third region of the stationary state for the 
parameter choice P > 2 , s > ^z^- Our aim is to derive a closed form of the time 
correlations of the slope at a fixed site in this region. The derivation is analogous 
to the one given in [5] for stochastically driven ordinary differential equations. 
Since in this region the slope is supercritical all the time the relaxation rule takes 
on a simple form. We simplify further by adjusting the left and right boundary 
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conditions. Together with the perturbation the simplified system can be regarded 
as a vector-valued stochastic process describable by a mathematically tractable 
stochastic linear recursion relation. Let 



i : No -^ R'^ 
(where ^(t) is the state at time t > 0) denote the process, and 

the initial state. Then the evolution equation becomes 

i(t) = Ai(t-i)+f(t) , t>o 

with A the relaxation matrix 



(24) 



(25) 



(2Q) 



I 



A 



1 



p 











\ 



p 



V 











1- ^ 

p 



p 



p 



1-^ 



{N xN - matrix) (27) 



? / 



and r(i) the stochastic perturbation. This relaxation differs from the one considered 
in section 3.2 through the boundary condition for the first and last site: The 
relaxation matrix here ensures conservation of the overall sum of slopes in the 
vector ^. The random vector r(i) has iV — 1 equally probable outcomes, 



f«:=(0,...,0,-s, ^ ,0,...,0) , ^ = 2,...,iV . 
i.th comp. 

For different times the perturbation vectors are independent random vectors. 
The matrix A has eigenvalues 



f28) 



Xi. 



2 / (k- l)7r\ 

1 - p ( 1 - cos ^ ^ ^ ) , A; = l,...,iV 



f29) 



2.(1) li.^) 

with corresponding eigenvectors v ,...,« , 



v'P 






1 



^=l,...,iV 



2 (fc-l),r(2»-l) 



A; = 2,...,iV , i = l,...,iV 



f30) 
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which we combine in the matrix 



,-Sj)\n 



W := (.«,,)-=, := (^rO" =1 



f311 



Note that Ai = 1 and |Afc| < 1 for A; > 1 , P > 2. 

In the following we want to derive out of the known statistical properties of the 
perturbation process T the statistical properties of the process ^. In particular, we 
will calculate the time correlation function and the spectral density function. The 
process T has the averages 



<(Tit))i> =0 for i = 2,...,iV- 1 



< (m)N > 



N-l 



f32) 



and the variances 



< [(r(ox- < (rm >] mtyh- < (m). >] > 

< imumh >-< imh >< (mh > 

, i=j = lori = j = N 



m) 



N-i (N-iy 

N-l 
N-l 



(Ar-1)2 





, K-il = 1 

, i = l^j = Nori = N,j=l 
, else 



(34) 



A central role in the derivation is played by the formal solution of the recursion 
relation f26): 



y t>0 : iit)=J2 ^'""f (m) + A'c . 



m) 



The powers of A can be expressed in the form 



N 



iA^),, = J2W,,W,,Xl . 



m) 



k = l 



Insertion of (36) in (35) and calculation of the expectation value leads to 

t N 



< 6(0 > 



N 



xEE^'^^i^^l" 



m=l k=2 
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t N 

m=l k=2 

N 

+ J2 W.kW,,Xlc, . (37) 

j,k=i 

Since we are interested in the stationary state we take the limit t ^ oo: 
< ^i(oo) > := lim < Ci{t) > 



Ps . PsN +1 1 ^ 



^c, . (38) 



N -I 2 N -I N . ^ 

A similar calculation can be done for the second central moments of the random 
vector S,{t): 

<y.,{i) :=< (6(0- < 6(0 >)(6-W- < 6-W >) > • (39) 

In the stationary state we find 

(Tij{oo) := ^liij|^o-ij(0 (40) 

= E E 9k,kMr,Wk,r,W,r,Wk,r, _ . . ■ (41) 

The temporal correlation functions of the process ^ are defined by 

K,,{T,t) : = < (6(i + r)- < at + r) >) {i,{t)- < ^.(t) >) > (42) 

for t > , T £ Z such that t + T > 0. In the stationary state they can be calculated 
to be 



Vr>0 : Kijir) := lim Kij(T,t) 

N N 

= J2J2^^rW,A>kAoo) (43) 

r = 2 fc = l 

VrGZ : A',,(r) = A',,(-r) (44) 

The spectral density function is the Fourier transform of the autocovariance 
function, 

= ii ('«(-) + 1 1 K>Hv..,(oo) ,!^^—; ^y . (45, 

The autocovariance function and spectral density function calculated from (43) 
and (45), respectively, are compared with simulation results in Fig. 9. The simu- 
lated data cannot be extended to very low frequencies because of time and memory 
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constraints on the time series analysis. The calculated spectral density function 
shows the two different regions of algebraic behaviour as mentioned in section 2.3. 
For very low frequencies it approaches the value Sjj(O). 

The mechanism for the short time correlations which correspond to the higher 
frequencies is the damping of perturbations. For P > 4 all eigenvalues of the 
relaxation matrix are positive. The one with the smallest value 

2 / (N - 1)tt\ 4 , , 

Xn = 1-j(i-cos ^ ^ ^ j ^ 1 - _ (46) 

corresponds to the fastest decay and may thus serve as a measure for the relaxation 
rate. An excitation of this eigenvalue will decay as (1— 4/P)*. The typical relaxation 
time T can then be estimated from (1 — 4:/P)'^ = 0.01 (decay to 1%) or T = 
ln0.01/ln(l — 4:/P). This time corresponds to the frequency 

oj = - = ln(l )^- ]=— (47) 

T InO.Ol ^ P^ InO.Ol \ Pj P ^ ^ 

which is quite a good estimate of the transition frequency between the two scaling 
regions. For instance for the parameter choice used in Fig. 9 this estimate predicts 
the second transition frequency to be near loguj = —0.5. 

We suspect perturbations at different time steps to be the mechanism respon- 
sible for the long time correlations which correspond to low frequencies. That the 
low frequency cut off in the spectral density function is due to a finite size effect 
can be seen through inspection of eq. (45): Sjj is a superposition of functions of 

the form 

„ , , 2A cosu) — 2A^ , -, r / X 

For small frequencies this is a Lorentzian 

, , , 2A(1- A)- Acj2 ^ ^ 

The Lorentzian in eq. (45) with the highest maximum and the narrowest width is 
the one corresponding to the eigenvalue closest to 1, i.e. 



2, 


vr X 


1 ^2 


-F^'- 


-cos— ) ^ 1- 




PiV2 



^2 = i-p(i--^]^) = i-p]^ • (50) 

The transition frequency uj^ may be estimated by that frequency where the domi- 
nant Lorentzian drops to a fraction c of its maximum, i.e. 

/a.(^c) = c/a,(0) , cG]0,l[ . (51) 

After a few steps one arrives at: 



_ / 2(l-A)2(l-c) 
^"V (l-A) + 2cA 
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(1-A)i 


J'-' 


1 vr^ 1 

PN'-'y 


ll-c 
c 



(52) 



The choice c = 0.95 leads to an estimate for the low frequency cut off in Fig. 9 
of loguj = —3.4. A plot of eq. 52 for increasing lattice size N shows a decrease in 
the estimated transition frequency which supports the assumption of a finite size 
effect . 

5 Summary and conclusions 

In this paper we have analyzed a stochastically driven running one- dimensional 
sandpile automaton by looking at the dynamics on a finite lattice. We could explain 
qualitatively the stationary slope profiles and calculate the transition sites between 
the different phases. The profiles are similar to the ones of the limited local sandpile 
model found in [3]. The limited local model lacks a linearly increasing part since 
it does not incorporate weighted relaxation. 

The kind of transport argument used here is well adapted to discrete time 
dynamics on a finite lattice and allows to characterize the stationary state without 
looking at a finite size scaling theory for the corresponding continuum model. It also 
helps to avoid problems met when approximating the dynamics by a continuous, 
nonlinear diffusion equation as encountered in [3]. In particular, the ability to 
achieve the required transport via undamped running waves seems to be the origin 
of the singularity of the diffusion constant in such models. 

We calculated exactly the nontrivial scaling autocovariance functions and spec- 
tral density functions of time series recorded in phase III using well known methods 
for stochastically perturbed linear recursion relations. The final form for the spec- 
tral density function, eq. (45), allows to estimate the transition frequencies between 
the scaling regions. It also contains the scaling behaviour but we have not succeeded 
in estimating exponents analytically. 

How much of these observations and methods generalizes to higher dimensional 
automata is not obvious. One might expect to find spherical waves instead of our 
linear waves, but as the wave propagates outwards the mass it transports is diluted 
up to the point where it stops propagating. This transition between running and 
stationary sites complicates the analysis. Similarly, because of the dilution it is 
difficult to build up the linear, always supercritical profile, in which one could 
apply methodes for stochastically driven linear maps. The easiest way to achieve 
this is an 'inverted sandpile', where sand is drained in the center and not at the 
circumference of the pile. Since the mass is then transported towards the center, 
the fiux rates increase and the sites will always be supercritical. This might serve 
as a model for the formation of a hydrodynamic drainage network (see for instance 
the experiment described in [6]). 
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Figure 1: Schematic picture of the linear automaton before (a) and after (b) a perturbation, hi denotes 
the height at site i and Cc the slopes. The values after a perturbation are marked by *. 

Figure 2: Division of the P — s— plane. The crosses mark the parameter values used in the simulation and 
are accompanied by references to figures which contain the stationary slope profiles. The two continous lines 
divide the plane into three different parts characterized by distinct structures of the profiles (see the text for 
details). They are predictions which are obtained by transport arguments described in subsection 3.3. The 
two dotted lines show the divergence points of the curves. For parameters chosen out of the shaded region 
the automaton does not converge to a stationary state. 

Figure 3: Expectation value of the slope profile for the parameter choice _P = 6 , s = 0.25 and snapshots of 
the time evolution. The profile is a nearly constant subcritical line. Subcritical wave packets are characteristics 
of the dynamics in all parts of the lattice. 

Figure 4: Expectation value of the slope profile for the parameter choice _P = 6 , s = 1.0 and snapshots 
of the time evolution. Compared to Fig. 3 one observes a second regime in the slope profile: A regime with 
supercritical slope in which the dynamic is dominated by supercritical wave packets. The two parts of the 
profile are nearly horizontal. 

Figure 5: Expectation value of the slope profile for the parameter choice _P = 6 , s = 2.5 and snapshots of 
the time evolution. Compared to Fig. 4 this profile shows more structure: There is an additional third part 
which increases linearly and in which perturbations are damped. 

Figure 6: Expectation value of the slope profile for the parameter choice P = 1.5 , s = 1.0 and snapshots 
of the time evolution. This pair of parameters lies in the same region of the phase diagram in Fig. 2 as the 
one used in Fig. 3. The only difference is a slight ascent in the last part of the lattice. 

Figure 7: Expectation value of the slope profile for the parameter choice P = 1.5 , s = 2.5 and snapshots 
of the time evolution. These parameters are taken out of region II of the phase diagram in Fig. 2; one has 
to compare this profile with the one given in Fig. 4. The dynamics observed in sub- and supercritical parts 
agrees, while the plateaus in Fig. 4 degrade to a monotonically increasing profile. 

Figure 8: Autocovariance function (a) and spectral density (b) for a time series recorded at site 160 in the 
automaton with parameters _P = 10 , s = 2. The site 160 lies in the linearly increasing part of the stationary 
profile. The two scaling regions in the spectral density have slopes -0.6 and -1.44, respectively. 
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Figure 9: Measured and calculated (from (43) and (45)) autocovariance function (a) and spectral density 
function (b) of the slope at site 25 in the model described by the recursion relation (26) with parameters 
_P=10, s = 2 , N = 50. The calculated curves are shifted by 0.5 units of the vertical axis. The two scaling 
regions in the calculated curve have scaling exponents of -0.53 and -1.45, respectively. 
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